// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

#ifndef _PMATHS_H_
#define _PMATHS_H_

/*! 
* \defgroup pmaths Pmaths Utility Library
* \brief	Maths utility library containing sub-routines obtained from PLINK.
* \details	The math utility library is lifted 'as is' from PLINK.<br>
*			Other than a few type conversions and changes to the file header file import declarations
*           to allow the code to compile with VC++, the algorithmic sections of the code remain unchanged.<br><br>
*           N/B Large sections of these codes seem to derive from Numerical Recipes 3 (NR3).
*/

/*!
 * \file
 * \ingroup pmaths
 * \brief	PMATHS Utility Library
 * \details Contains sub-routines lifted from PLINK that are associated with the hard-sums
 *			part of the statistical calculations.
 */
#include "types.h"

#ifdef __cplusplus
extern "C" {
#endif

/*! 
\brief Maths Return enumeration
\deprecated Mostly ununsed, will be removed after Epistasis (QT) refactored.
*/
enum MATHS_RETURN {MATHS_FAILURE = 0, MATHS_SUCCESS = 1, MEMORY_FAILURE = 2};

/*!
\brief check a number for any divide by zero errors.
\return TRUE if value okay
*/
BOOL realnum(double d);

/*!
\brief Cumulative Distribution Function CHI-Square distribution wrapper
\param[in] x the variable
\param[in] df Degree of Freedom
\return double or -9 on failure
*/
double chiprobP(double x, double df);

/*!
\brief Get list of randomised array indices.
\details
	The array is filled with values 0 to nA-1 and
	the order then modified by calls to get_rand();
\param[in, out] a Source array
\param[in] nA Size of the source array
\return 1 if successful, 0 on failure
\see get_rand()
*/
int permute(int *a, int nA);

/*!
\brief Get a random number
\details Get a random integer from 0 to n-1.
\param n Maximum random value
\return integer from 0 to n-1
\see genrand_real2()
*/
int get_rand(int n);

/*!
\brief The Logarithm Of The Gamma Function
\details Appears to be from NR3
\param xx Value to pass to the distribution function
\return double
*/
double gammln(double xx);

/*! 
\brief Lower tail quantile for standard normal distribution.
\param p Value to pass to the distribution function
\return double
*/
double ltqnorm(double p);

/*! \brief Evaluates the continued fraction for the incomplete beta function. */
double betacf(const double a, const double b, const double x);

/*! \brief The Incomplete Beta Function */
double betai(const double a, const double b, const double x);

/*! \brief Unknown */
double pT(double T, double df);

/*! \brief Calculate Pythagoras distance between 2 vectors. */
double pythag(const double a, const double b);

/*! \brief Hardy-Weinbury Equilibrium
\details
	This function implements an exact SNP test of Hardy-Weinberg
	Equilibrium as described in Wigginton, JE, Cutler, DJ, and
	Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg
	Equilibrium. American Journal of Human Genetics. 76: 000 - 000
\author Jan Wigginton */
double SNPHWE(int obs_hets, int obs_hom1, int obs_hom2);

/*! \brief Normal Distribution */
double normdist(double z);

/*! \brief Co-linearity Check 
\details	PLINK utility function concerned with performing co-linearity checks of a 
			matrix associated with a linear algebra calculation.
*/
BOOL check_vif(MATRIX_T X, int x_rows, int x_cols);

/*! \brief Perform Single-Vector decomposition on a input matrix */
int svd_inverse(MATRIX_T u, MATRIX_T r, int nRows, int nCols);

/*! \brief SVD Compare (NR2) */
int svdcmp(MATRIX_T a, int a_rows, int a_cols, VECTOR_T w, MATRIX_T v);

/*! \brief Matrix multiplication 
\details c = a * b
\param a Matrix a
\param ar No. rows of a
\param ac No. columns of a
\param b Matrix b
\param br No. rows of b
\param bc No. columns of b
\param c Results of matrix multiplication
\return 1 on succees
*/
int mult_matrix(MATRIX_T a, int ar, int ac, MATRIX_T b, int br, int bc, MATRIX_T c);

/*! \brief Perform SVD Back-substution */
enum MATHS_RETURN svbksb(MATRIX_T u, int u_rows, int u_cols, VECTOR_T w, MATRIX_T v, VECTOR_T b, VECTOR_T x);

/*! \brief For integer representing array of size n, get number of unique 'pairwise' interactions */
long getPairwiseTests(long n);

/*! \brief For a given number of elements in an array, calculates the number of pairwise tests to run on each node. */
long getPairwiseTestsPerNode(long nElements, long nNodes);

/*! \brief Calculate the epistasis indices for the non-permuted calculation. 
\details 
	For a input of array of size N (probably SNPs), the cluster of size nodes and the 
	the processor MPI rank; this function calculates the outer loop limits for the 
	epistasis algorithm so that the non-permuted calculation can be distributed across
	N-cores.
\param[in] N Size of an input SNP array
\param[in] nodes No. cores allocated to a MPI calculation
\param[in] mpi_rank The local MPI rank of a processing thread
\param[in,out] indices Array where the calculated indices are stored.
\author David Hoyle
\return TRUE if successful
*/
BOOL getDHLoopLimits(int N, int nodes, int mpi_rank, int *indices);

/*! \brief Get the P-value associated with an CHISQ value. */
double getChisqPvalue(double chisq);

/*! \brief Update permutation counter arrays for 1-D CHISQ calcuation.
\param[in] size Size of the input array(s).
\param[in] results test statistic for the permuted calculation.
\param[in] originals test statistic for the non-permuted calculation.
\param[in,out] R Counter for the empirical p-value calculation (EMP1)
\param[in,out] maxR Counter for the corrected family-wise p-value calculation (EMP2)
\return 1 if successful
*/
int perm_update(int size, double *results, double *originals, int *R, int *maxR);

/*! \brief Randomises a sequential integer array.
\details 
	Randomises a sequential integer array. Unlike \ref permute() this function
	does not populate the source array with the ordered indices. 
	This means that the array is constantly 2 re-shuffe at each iteration iteration.
\param[in,out] src Pre-filled source array of integers
\param[in] nSrc Size of the source array
\see permute()
\see fill_indices()
*/
void shuffle(int *src, int nSrc);

/*! 
* \brief Get the number of tasks allocated to a specific number of cores.
* \details Tasks are either a slice of data from an input array or the number of 
*          permutations that must be run as part of a statistical calculation.
* \param[in] numprocs Number of processors/cores allocated to a job.
* \param[in] nSample Total number of tasks allocated to a job which are to be sub-divided 
*            amongst the processors.
* \return This method returns a 'floored' unsigned integer
* \see _get_tasks_per_processor()
*/
unsigned long get_tasks_per_processor(unsigned int numprocs, unsigned long nSample);

/*! 
* \brief Get the number of tasks allocated to a specific number of cores.
* \details Tasks are either a slice of data from an input array or the number of 
*          permutations that must be run as part of a statistical calculation.
* \param[in] numprocs Number of processors/cores allocated to a job.
* \param[in] nSample Total number of tasks allocated to a job which are to be sub-divided 
*            amongst the processors.
* \return This method returns a 'ceiling' unsigned integer for even section lengths required by MPI_Gather().
* \see get_tasks_per_processor()
*/
int _get_tasks_per_processor(int numprocs, int nSample);

/*!
\brief Calculate a odds ratio.
\details
<p>Method provided by Wolf (<em>see <a href=http://slack.ser.man.ac.uk/theory/association_odds.html>Odds Association</a></em>).<br>
Based on asymptotic, large %sample approximation.</p>
\param [in] A1 Affected count one
\param [in] A2 Affected count two
\param [in] U1 Unaffected count one
\param [in] U2 Unaffected count two
\param [in,out] OR Pointer to store the Odds ratio
\param [in,out] lowerCI Pointer to store the lower confidence interval
\param [in,out] upperCI Pointer to store the upper confidence interval
\param [in] zt Normalised confidence interval threshold (zt = ltqnorm(1 - (1 - ops.ci_level) / 2))
\returns 1 or success, 0 on failure
*/
int getOddsRatio(int A1, int A2, int U1, int U2, double *OR, double *lowerCI, double *upperCI, double zt);

#ifdef __cplusplus
}
#endif

#endif // _PMATHS_H_

